/*
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

// -- matriz.h
double** alocar_matriz(int altura, int largura);
void desalocar_matriz(double **M, int altura);
double** ler_matriz(FILE *file, int *J_ORDER, int *J_ROW_TEST, double *J_ERROR, int *J_ITE_MAX);
double* ler_b(FILE *file, int *J_ORDER);

// -- main.c
int main(int argc, char **argv)
{
    int J_ORDER = 0, J_ROW_TEST = 0, J_ITE_MAX = 0;
    double **A, *B, *D, *aux;
    double J_ERROR, ctrl, error_ctrl = 1, soma = 0, test = 0, max = 0, max_diff = 0;
    double dtest = 0;
    int i, j, k;
    FILE *file;

    if(argc != 2){
        printf("usage: jacobiSeq <filename>");
        exit(-1);
    }

    file = fopen(argv[1], "r");
    if(file == NULL){
	printf("Erro ao abrir o arquivo\n");
	exit(-1);
    }

    A = ler_matriz(file, &J_ORDER, &J_ROW_TEST, &J_ERROR, &J_ITE_MAX);             //Matriz principal
    B = ler_b(file, &J_ORDER);                                //Vetor B
    D = (double*) calloc(J_ORDER , sizeof(double));     //D vetor com os valores encontrados na iteracao atual (inicialmente 0,0,0...0)
    aux = (double*) calloc(J_ORDER , sizeof(double));   //auxiliar para o D da proxima iteracao

    for (k = 0; k < J_ORDER; k++) // "divide a matriz pela diagonal" e zera a diagonal
    {
    	ctrl = A[k][k];
    	if(ctrl == 0){
            printf("Matriz invalida");
            exit(-1);
            }
    	B[k] = B[k]/ctrl;        //divisao dos componentes de B pelo respectivo elemento da diagonal (B torna-se B*)

    	for (i = 0; i <J_ORDER; i++)
			A[k][i] = A[k][i]/ctrl;   //divisao dos componentes da linha de A pelo respectivo elemento da diagonal (A torna-se A*)

        A[k][k] = 0;   //zera a diagonal, mais eficiente do que uma comparacao dentro do loop
    }

    for (k = 0; k < J_ITE_MAX && error_ctrl; k++) //enquanto k < numero de iteracoes maximo e erro nao satisfeito
    {
        max = 0;
		for (i = 0; i < J_ORDER; i++){     //algoritmo de jacobi
            soma = 0;
            test = 0;
            for (j = 0; j < J_ORDER; j++){
                soma += A[i][j] * D[j];
                test += fabs(A[i][j]);
            }

            aux[i] = B[i] - soma;  //uso de vetor auxiliar para nao alterar o D enquanto ainda esta sendo usado
            if (test >= 1){
               printf("Sistema divergente\n");
               return -1;
            }
        }
        max_diff = 0;
        for (i = 0; i < J_ORDER; i++){

            if(fabs(aux[i]) > max)
                    max = fabs(aux[i]);  //max pega o maior entre os valores de aux[i] (valores parciais do vetor resposta B atualmente)
            if(fabs(aux[i]) - fabs(D[i]) > max_diff)
                max_diff = fabs(aux[i]) - fabs(D[i]);  //pega a maior das diferencas entre cada aux[i] e D[i]

            D[i] = aux[i]; //copia o auxiliar para D para atualiza-lo para a proxima iteracao de k
        }
        if(max != 0 && max_diff / max < J_ERROR) //checa se satisfaz o erro
            error_ctrl = 0;  //se sim zera error_ctrl que quebra o loop
    }


    for (i = 0; i < J_ORDER; i++){ // volta a diagonal para 1 (estava com 0 anteriormente para o algoritmo de jacobi)
    	A[i][i] = 1;
        dtest += D[i] * A[J_ROW_TEST][i];
    }

    printf("\n---------------------------------------------------------\n");
    printf("Iterations: %d\n", k);
    printf("RowTest:%d => %lf =? %lf\n", J_ROW_TEST, dtest, B[J_ROW_TEST]);
    printf("---------------------------------------------------------\n");

    //Limpando
    desalocar_matriz(A, J_ORDER);
    if(B != NULL)
        free(B);
    if(D != NULL)
        free(D);
    if(aux != NULL)
        free(aux);
    fclose(file);

    return 0;
}

//Aloca matrizes
double** alocar_matriz(int altura, int largura)
{
    int i;
    double** M = (double**) malloc(altura * sizeof(double*));
    for(i = 0; i < altura; i++){
        M[i] = (double*) malloc(largura * sizeof(double));
        if(M[i] == NULL){
            printf("Erro na alocacao da matriz\n");
            exit(-1);
        }
    }
    if (M == NULL)
        exit(-1);
    return M;
}

//Desaloca matrizes
void desalocar_matriz(double **M, int altura)
{
    int i;
    for(i = 0; i < altura; i++)
        if(M[i] != NULL)
            free(M[i]);
    if (M != NULL)
        free(M);
}

//Leitura e alocação de matrizes (chama alocar_matriz)
double** ler_matriz(FILE *file, int *J_ORDER, int *J_ROW_TEST, double *J_ERROR, int *J_ITE_MAX)
{
    fscanf(file, "%d\n", J_ORDER);
    fscanf(file, "%d\n", J_ROW_TEST);
    fscanf(file, "%lf\n", J_ERROR);
    fscanf(file, "%d\n", J_ITE_MAX);
    double **M = alocar_matriz(*J_ORDER, *J_ORDER);
    int i, j;

    for (i = 0; i < *J_ORDER; i++)
        for (j = 0; j < *J_ORDER; j++)
            fscanf(file, "%lf", &M[i][j]);

    return M;
}


//Leitura do vetor B
double* ler_b(FILE *file, int *J_ORDER)
{
    int i;
    double *B = (double*) malloc (*J_ORDER * sizeof(double));
    for (i = 0; i < *J_ORDER; i++)
        fscanf(file, "%lf\n", &B[i]);
    return B;
}
